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Abstract 

We calculate the second order viscous correction to the kinetic distribution, 5/(2) , and use this 
result in a hydrodynamic simulation of heavy ion collisions to determine the complete second 
order correction to the harmonic spectrum, v n . At leading order in a conformal fluid, the first 
viscous correction is determined by one scalar function, xop- One moment of this scalar function 
is constrained by the shear viscosity. At second order in a conformal fluid, we find that 5f(p) can 
be characterized two scalar functions of momentum, xip an d X2p- The momentum dependence 
Qh of these functions is largely determined by the kinematics of the streaming operator. Again, one 

"^ moment of these functions is constrained by the parameters of second order hydrodynamics, t w 

and Ai. The effect of 5/(2) on the integrated flow is small (up to V4), but is quite important for the 
higher harmonics at modestly-large pt- Generally, 5/(2) increases the value of v n at a given p?, 

r Z^ ] and is most important in small systems. 
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I. INTRODUCTION 

Motivated by the wealth of data on collective flow in heavy ion collisions, hydrodynamic 
simulations of these ultra-relativistic events have steadily improved. Today, event by event 
viscous hydrodynamic simulations reproduce the elliptic fow, the triangular flow, the higher 
harmonic flows, and the correlations between the flow harmonics and event plane angles 
[1]. The overall agreement with the data on this rich variety of observables constrains the 
shear viscosity of the QGP close to transition region. Because of this success, there is a 
current need to precisely quantify the systematic uncertainties in these simulations and the 
corresponding uncertainty in the extracted shear viscosity of the QGP [2]. This paper will 
address this need by computing the second order corrections to the thermal distribution 
function [3, 4] and by using these results to simulate the observed harmonic spectrum. 

Since the nucleus is not vastly larger than the mean free path, an important advance 
in hydrodynamic simulations was the inclusion of viscous corrections to the hydrodynamic 
equations of motion at first and second order in the gradient expansion. In particular Baier, 
Romatschke, Son, Starinets, Stephanov (BRSSS) determined the possible tensor forms that 
arise in the constituent relation of a conformal fluid at second order [3] . These equations 
are currently used in practically all viscous simulations of heavy ion collisions. It is satis- 
fying that the gradient expansion, which underlies the hydrodynamic approach, seems to 
converge [5], and seems to converge to the measured flow [1]. However, all simulations of 
heavy ion collisions make additional kinetic assumptions about the fluid at freezeout when 
computing the particle spectra [6]. Indeed, the phase space distribution of a viscous fluid 
f p (X) is modified from its equilibrium form, n p (X), by corrections at first and second order, 
$fp — ^/(i) + <5/(2)- Currently, all simulations of heavy ion collisions compute the viscous 
corrections to the distribution function at first order, while using second order corrections to 
the hydrodynamic equations of motion. The goal of this work is to remedy this inconsistency 
by computing viscous corrections to the distribution function to second order. Then, this 
second order correction is used to compute the harmonic spectrum, f n (pr)- As expected [7], 
these corrections are modest for small harmonic numbers, n, and small px, but grow both 
with n and p?- 

To compute the viscous corrections to the phase space distribution we will analyze kinetic 
theory of a conformal gas close to equilibrium in a relaxation time approximation. This 
extreme idealization is still useful for several reasons. First, a similar equilibrium analysis of 
QCD kinetic theory was used to determine the second order transport coefficients to leading 
order in a s [4]. This analysis (which will be discussed further below) makes clear that the 
details of the collision integral hardly matter in determining the second order transport 
coefficients. Indeed, we will see that the structure of the second order viscous correction 
5/(2) is largely determined by the kinematics of free streaming, rather than the details of the 
collision integral. Thus, an analysis based on a relaxation time approximation is an easy way 
to reliably estimate the size of such second order corrections in heavy ion collisions. Second, 
the overall normalization of 5/(2) is constrained by the second order transport coefficients 
r n and Ai, in much the same way that the shear viscosity constrains the normalization of 
5/(i). These constraints allow us to make a good estimate of the form of viscous corrections 
at second order for a realistic non-conformal fluid. 

With this functional form we study the effect of 5/(2) on the harmonic spectrum in heavy 
ion collisions. In Section II we outline how 5/(2) is computed in kinetic theory. Then, in 



Section III we discuss the practical implementation of this formula in a hydrodynamic code 
used to simulate heavy ion collisions. These results are used to simulate the linear response 
to a given deformation, e n . The linear response is a largely responsible for determining the v n 
in central collisions. In non-central collisions the linear response and the quadratic response 
determine the harmonic flow and its correlations [1, 8-11], but the quadratic response will 
not be discussed in this initial study. Finally, we will summarize the effects of <5/(2) m 
Section IV. 



II. 2ND ORDER CORRECTIONS TO THE PHASE SPACE DISTRIBUTION 

A. Notation 

Throughout the paper we will work with the metric rf" = (— , +, +, +), and d — 4 notates 
the number of space-time dimensions. Capital letters P,X denote four vectors, while lower 
case letters p, x denote three vectors in the rest frame. 

We are expanding around a fluid in equilibrium with energy density, pressure, and flow 
velocity equal to e(X), V(X), and U fl (X) respectively where U P U P = —1. Thus, the rest 
frame projector is A pu = g pu +U p U v , and the spatial derivatives and temporal derivatives are 
V M = A pv d v and D = u p d pi respectively. Finally bracketed tensors are rendered symmetric- 
traceless and orthogonal to w M , 

(A 1 *) = l - A^A 1 ; (A pa + A ap ) - — ^— A"" A pa A p(7 . (2. 1) 

Such tensors transform irreducibly under rotation in the local rest frame. A more elaborate 
example using a pu = 2 (V ,j,u v ) which appears in the algebra below is 



i°7«.A*2 <T M3A*4Jsym \ < fylM2 <T A*3/M/ + j . o t^Vl/M \ CT /*3 CT / U 4A/ J S y m 

"*" (d—l)(d+ 1) WM2 ^W/ sym ° ' V-Z) 

where the symmetrized spatial tensor is denoted with curly brackets: 

l °VlM2 ^MS^J sym = O rVlM2°7*3M4 + ° HlfJ>3° M3/*4 + 0/411*4^2^3 J ' K^'^J 

The equilibrium distribution function is n p = n(—P ■ U(X)/T(X)) where n(z) = l/(e z =F 
1), and f p {X) denotes the full non-equilibrium distribution. The rest frame integrals are 
abbreviated J = J d p/(27r) d-1 and primes (such as n') denote derivatives with respect 
to —P ■ U/T, so that n' = — n p (l ± n p ). The energy and squared three momentum in the 
rest frame are, E p = —P ■ U and p 2 = P p P u A flu , respectively. 



B. Hydrodynamics 

In evaluating Sf^) to second order we will need the hydrodynamic equations of motion. 
In the Landau frame these equations read 

T pu = eu p u u + PA"" + tc pu , T"" = , (2.4) 



where ir^ v is the viscous correction to the stress tensor. Throughout this analysis we are 
working with a conformal fluid, and consequently the bulk viscosity is zero ( = 0. For a 
conformal fluid the possible tensor forms of the gradient expansion for n^ through second 
order were established by BRSSS 
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+ Ai (a\a" x ) + A 2 (a\Q vX ) + A 3 (n\n vX ) + ..., (2.5) 
where —rja^ = —2r] (V^u u ) is the first order term, and the vorticity tensor is 

1 



W = A" a A^ (d a u? - d p u a ) . 



(2.6) 



The ellipses in Eq. (2.5) denote terms third order in the gradients. Using these equations of 
motion, the time derivatives of the energy density and flow velocity can be determined from 
the spatial gradients of these fields 
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(2.7a) 
(2.7b) 
(2.7c) 



In passing from Eq. (2.7b) to Eq. (2.7c) we have used the first order expression for 7r Ml/ = 
—rja^, the conformal temperature dependence of r\ oc T d_1 , and the lowest order equations 
of motion. 

In hydrodynamic simulations of heavy ion collisions the static form of the constituent 
relation Eq. (2.5) is not used. Rather, this equation is rewritten as a dynamical equation 
for 7r M j, which is evolved numerically [3] 
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(2.8) 



Similarly, when constructing Sf at first and second order we will systematically replace a^ u 
with —Tf^/r]. For fluids with an underlying kinetic description, the transport coefficients 
are additionally constrained, with A3 = and A2 = —2r]T n [4], and these relations used in 
the simulation. The appropriate values for Ai and r]T n will be discussed below. 



Kinetics 



To determine the viscous corrections to the distribution function we will solve the kinetic 
equations in a relaxation time approximation through second order in the gradient expansion 



f P = n p + Sf p = n p + (5/ ( i) + 5f {2) + . . 



(2.9) 



In a relaxation time approximation the Boltzmann equation reads 

P^f p (X) = ~ [f p (X) - n(-P ■ U.(X)/T.{X))] , (2.10) 

Up 

where the dimensionless coefficient C p is related to the canonically defined momentum de- 
pendent relaxation time 

C p = T 2 ^^-. (2.11) 

E p 

Following Ref. [7], we will parametrize the momentum dependence of the relaxation time as 
a simple power law 

T R <xE*- a , (2.12) 

with a between zero and one. As we will see below, a = gives a first order a viscous correc- 
tion which grows quadratically with momentum (which is known as the quadratic ansatz), 
while a = 1 gives a first order viscous correction which grows linearly with momentum (and 
is known as the linear ansatz). 

At leading order, the parameters T* and Ug which appear in the kinetic equation are 
equal to the Landau matched temperature and flow velocity, T and £/ M . However, starting 
at second order T* and U£ will differ from T and U^ by squares of gradients: 

T*(X) =T(X) + S%(X) , (2.13a) 

U?(X) =U^{X) + 6U?(X) . (2.13b) 

5T* and 8U£ will be determined at each order by the Landau matching condition: 

T^Uu = eu" . 

Expanding n(—P-U*/T*) we have 



n* p = n p + 5n* p 5n* p = n' p 



P-5U* 5% 



+ ..., (2.14) 



where we have used an obvious notation, n* = n(—P ■ U*/T*) and n p = n(—P ■ U/T). 

Then to determine Sf we substitute the expansion (Eq. (2.9)) into the relaxation time 
equation and equate orders. In doing so we use the hydrodynamic equations of motion 
through second order to write time derivatives of T(X) and C/ M (X) in terms of spatial 
gradients of these fields. For instance, using the equations of motion Eq. (2.7) and the 
thermodynamic identities 

de 1 _ J 1 \ „ 2 e + V 



and the relations p 2 = E 2 and c 2 = l/(d — 1) of a conformal gas, we have 
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E p P^A, lll2 d X 7c x ^ + \E 2 v c 2 sW 2 +.... (2.16) 



The term linearly proportional to a^ v is ultimately responsible for the shear viscosity, while 
the nonlinear terms contribute to 5/(2). 

With this discussion, we find that the Sf is determined by the hierarchy of equations: 

<% = C pK^^ > ( 2 - 17 ) 



and 
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p ' (e + p)T3 



E p P^ A^ 2 d X 7r^ + lE 2 p c 2 sV a 2 ] - ^P^Sffo. (2.18) 

Straightforward algebra uses the equations of motion to decompose Sfw into irreducible 
tensors, and determines the final form of SfZs and Sf7 2 ) (see Appendix A). 

We have put a superscript a in fZ, and 8f? 2 ) t° indicate that that we are using a^ u rather 
than n^ u in these equations. In realistic hydrodynamic simulations of heavy ion collisions 
71-^ is treated as a dynamic variable, and —r]a^ u is systematically replaced by n^. This 
yields the following reparameterization of Sf 



Sf(i) = - \C p n' p —7f^-^nv , (2.19a) 



, P^P" 

rjT 3 

piipv 

5/(2) =<J/(2) + \C P n p — ^- [7i> + 770^] , (2.19b) 



where we have replaced <t m „ with —Ti^/rj in the first order result, and appended the difference 
between these two tensors to the second order result so that, 5/(i) + 5/(2) = fifm + 5/£) U P 
to third order terms. 

To record the result for 5/(2), we first review the familiar first order case. At first order, 
5/(i) is described by a dimensionless scalar function xop(E p /T) 

pinpm 
5/i = Xo P — j^- 7W 2 , (2.20) 

which has been extensively studied in the literature, and determines the shear viscosity [7]. 
In the relaxation time approximation this function is related to the relaxation time 

Xv P = -\C p n' p . (2.21) 

One moment of this function is constrained by the shear viscosity. Indeed, from the defining 

relation 
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Sf P , (2.22) 



we determine the shear viscosity 



* = id-i)(d + i)T'li x »- (2 - 23) 

and a constraint on 5/( 2 ) 

piipv 

v 



This constraint reflects the reparameterization of <j^ v in the first order 8f(r\ with -k^ '. For 
later use and comparison, we note that the enthalpy is 

which can be obtained by comparing the stress tensor from kinetic theory for small fluid 
velocities (i.e. u^ ~ (1, v) with v < 1) to ideal hydrodynamics, T° l ~ (e + V)v l [12]. 

At second order the function 5f^) is described by two dimensionless scalar functions Xip 

and X2 P 

Xi P = - \C P Xop , (2-26) 

X2 P =C P Xop • (2.27) 

Two moments of these scalar functions are constrained by the second order transport coef- 
ficients r]T n and Ai 

Al + ^ =(d-l)(d + 8 l)(d + 3)Tt I ** Y p ' (2 ' 28) 



(d-l)(d+l)T 5 
Then the functional form of 5f(2) is 
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(2.29) 



Yl P^l P^2 pM3 p/^4 y, pilipH2pH3 

5f( 2 ) =-f y& (VAw) + -T- f~ 5 i( d + 2) (^i W V w InT) - (V M1 7r M2 ^ 3 )] 

+ ^ 2 > (2-30) 

where the four scalar functions £i p , ^2 P , ^3 P , £ip are linearly related to xop, Xip 5 X2 P 

gip = Xip / , , Q N ^^(^tt + Ai), (2.31a) 

(d + 3) ^ 

e 2p = |^^ P -Xo P , (2.31b) 

1 7V 



?3p = -X2 P ^ , n + 2xop7^ P + a P ,n p , (2.31c 
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(d+1) + xo v 

<Up = Xip ^_ 1 w rf 1 n - ^2p7j— pr - Xop--^ p c s + a Et n p ti p , (2.31c!) 



with p = p/T and _E p = E p /T. 

The coefficients oe* and a p „, come from Eq. (2.14) and are adjusted so that the Landau 
matching conditions are satisfied. More specifically, we choose 5U* and 8T* in Eq. (2.14) so 



that 

P ■ 5U P Ml 
Y^ = a ?* ^3 [A w a 2 9 Ai tt AiA2 ] , (2.32a) 

-E v ^=a Et E v ^- 2 . (2.32b) 

Then integrating over f p (X) to determine the stress tensor and demanding that Eq. (2.24) 
(which is a restatement of the Landau matching condition), we conclude that 



1 + d) (-) , (2.33a) 
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a °-='-^-j^'-ir + m^)[i) ■ < 2 - 33b > 

Despite being somewhat complicated, the functional form of 5/(2) * s severely constrained, 
and is bounded by the transport coefficients r/, Ai and r/r^ through Eqs. (2.23), (2.28), and 
(2.29). For a single component classical gas with the quadratic ansatz a = 0, Eq. (2.23) 
shows that C p = - , and the three scalar functions which determine 5/(2) can be simplified to 

77 1 fTI\ 2 1 /"l]\ 2 

Xo P = 7T s n P ' Xip = i {-) n v , X2 P =7> {-) n P ■ (2-34) 

We will discuss the implementation of (5/(2) in the next section. 



III. IMPLEMENTATION IN SIMULATIONS OF HEAVY ION COLLISIONS 

In this section we will implement the 5/(2) corrections in a 2+1 boost invariant hydrody- 
namic code. A full event-by-event simulation of heavy ion collisions with 5/(2), together with 
a comparison to data, goes beyond the scope of this initial study. Nevertheless, the effect 
of 5/(2) in larger simulations can be anticipated by understanding how the linear response 
is modified by 5/(2). Indeed, the qualitative features of event- by-event hydrodynamic simu- 
lations of heavy ion collisions (including the correlations between the harmonics of different 
order) are reproduced by linear and quadratic response [8-10]. In central collisions the linear 
response is sufficient, and was recently used to produce one of the best estimates of the shear 
viscosity and its uncertainty to date [2]. We will calculate the linear response to a given 
deformation e n in order to estimate the influence of 5/(2) on v n . 

In a given heavy ion event, the particle spectrum is expanded in harmonics 

f/A ' (/A ' v ^i; n (p r )e in( ^-*" (pr)) 4 ( niuulrx (cnj. I . (3.1) 
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where <p p is the azimuthal angle around the beam pipe. In a linear response approximation 
the n-th harmonic, v h (pt), in the event is assumed to be proportional to the n-th cumulant, 
e n , which characterizes the deformation of the entropy distribution 1 . Specifically, we first 
define the normalized entropy distribution at an time initial time tq 

_ t s(x ± ) 

P{Xl_) = y-^ r • (3.2) 



1 Note that we will use cumulants rather than moments to characterize the deformation [10, 13]. 



Writing the coordinates in the transverse plane x±_ = (x, y) as a complex number, z = 
x + iy = re 1 '*', we define the first six cumulants characterizing the harmonic deformations of 
the initial distribution 
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(3.3a) 
(3.3b) 
(3.3c) 
(3.3d) 
(3.3e) 
(3.3f) 



where (. . .) denote an average over p(x±). In a linear response approximation the orientation 
angle of n-th harmonic ^nipr) is aligned the cumulant angle $ n . 

The linear response coefficient v n (jp>r)/e n is independent of many of the details of the 
initial state [13], and can be reasonably computed by initializing 2 + 1 boost invariant 
hydrodynamics with a deformed Gaussian distribution, where the rms radius and amplitude 
of the Gaussian are adjusted to match the rms radius and total entropy of the event. For 
example, to simulate v s at RHIC at an impact parameter of b = 7.45 fm we initialize a 
Gaussian deformed by 63 = 0.05 with $3 = 



t s(x,t ) = C s {N p ) 



1 + 



(r 3 ) e 3 
24 





r 2 



7r (r 2 ) 



(3.4) 



where C s sets the total multiplicity in the event. Here (r 2 ) and (N p ) are computed using 
the Phobos Glauber model [14]. In order that the total entropy closely matches the total 
entropy in more complete simulations [5, 15], we take C s = 15.9 and 28.04 at RHIC and 
LHC respectively. This procedure has been used previously by the authors to determine the 
linear and non-linear response [10, 11, 13]. 

After initializing the Gaussian, we evolve the system with second order hydrodynamics, 
Eq. (2.4) and Eq. (2.8), using a variant of the central scheme developed previously [16, 17]. 
Then for a given n-th order harmonic perturbation e n we compute v n /e n by performing the 
freezeout integral at a constant temperature. This evolution requires an equation of state 
and specified hydrodynamic parameters at first and second order. In what follows we will 
consider a conformal equation of state for a single component classical gas, and a lattice 
motivated equation of state previously used by Romatschke and Luzum [5] . 

Since it is only for the conformal equation of state p = |e that the analysis of Section II 
is strictly valid we will discuss this case first, and then discuss the necessary modifications 
for a lattice based equation of state. To keep the final freezeout volume of the conformal 
equation of state approximately equal to the much more realistic lattice based equation of 



Momentum Dependence 

Linear Ansatz (a = 1) 

Quadratic Ansatz (a = 0) 



nr^ 
(d+l) = 5 

(d + 2) = 6 



Ai 

(d+l) 2 /(d + 3)=25/7 

(d + 2) = 6 



TABLE I. A compilation of rescaled second order transport coefficients for a linear and quadratic 
ansatz in a relaxation time approximation for classical statistics [4]. All numbers in this table 
should be multiplied by n 2 /(e + V). In a relaxation time approximation A2 = —2nT n and A3 = 
[4]- 



state, we choose the final freezeout temperature (T fo = 96 MeV) so that the entropy density 
at freezeout Sf rz = 1.87 fm~ 3 equals the entropy density of a hadron resonance gas at a 
temperature of T = 150 MeV. The relation between the temperature and energy density 
for the conformal equation of state is e/T 4 = 12.2, which is the value for a two flavor ideal 
quark-gluon plasma. The motivation for these choices, the parameters of the conformal 
equation of state, and further details about the initial conditions and freezeout we refer to 
our previous work - see especially Appendix B of Ref. [13]. 

The second order transport coefficients rjT n and Ai are all constrained by the momentum 
dependence of the relaxation time and the shear viscosity through Eqs. (2.23), (2.28), (2.29). 
As discussed in Ref. [7] , there are two limits for this momentum dependence which span the 
gamut of reasonable possibilities. In the first limit the relaxation time grows linearly with 
momentum, and a = in Eq. (2.12). This is known as the quadratic ansatz, and is most 
often used to simulate heavy ion collisions. In a similarly extreme limit the relaxation 
time is independent of momentum, and a = 1 in Eq. (2.12). This is known as the linear 
ansatz, and this limit provides a useful foil to the more commonly adopted quadratic ansatz. 
Once the shear viscosity and the momentum dependence of the relaxation time are given, 
the collision kernel is completely specified in the relaxation time approximation, and all 
transport coefficients are fixed. For a linear and quadratic ansatz we record the appropriate 
second order transport coefficients in Table I. 

So far this section has detailed the initial and freezeout conditions, as well as the second 
order parameters which are used in the conformal equation of state. Fig. 1 shows the re- 
sulting elliptic flow for the quadratic and linear ansatze for a conformal equation of state 
including the first and second order 5f. A conformal equation of state has a strong expan- 
sion, and, as a result, generally over estimates the magnitude of the 5f correction. Thus 
the conformal analysis provides a schematic upper bound on the magnitude of the 5/(2) 
correction. Further discussion of these results is reserved for Section IV 

Strictly speaking the analysis of Section II is useful only for a single component conformal 
gas. Nevertheless, we believe the usefulness of the analysis extends beyond this limited 
regime [4]. Indeed, examining the steps in the derivation one finds that only very-few non 
conformal terms appear at each order. For instance, if non-conformal corrections are kept 
in Eq. (2.17) one finds 
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FIG. 1. Differential V2(pt)/^2 f° r RHIC and LHC initial conditions, for a linear and quadratic 
ansatz, and a conformal equation equation of state (CEOS). The n p curves shows the flow from 
second order hydrodynamics without the viscous correction to the distribution function; Sfn) and 
(5/(!) + <5/(2) show the flow with the viscous correction at first and second order; and finally, the 
5/m result uses — r/o" MI/ instead of tt^ in the first order result (see Eq. (2.17)). The freezeout 
temperature is chosen so that the freezeout entropy density of the conformal gas equals that of a 
hadronic resonance gas at a temperature of T = 150 MeV. 



which shows that non-conformal terms (the second term in Eq. (3.5)) are either suppressed 
by | — c^, or are suppressed at high momentum relative to the conformal terms. 

To extend our analysis to a multi-component non-conformal equation of state we have 
followed the simplified treatment that is used in almost all simulations of heavy ion collisions. 
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First, we will treat all species independently 

p»d,r P {x) = ~ [r P {x) 



n 



■P-U m (X)/T,(X))] 



(3.6) 



where a = it, K, p, . . . is a species label 2 . We will also adopt the quadratic ansatz a = 0, so 
that Cp is independent of momentum. Then for every species we define the partial entropy 
and shear viscosity as in Eqs. 2.23 and 2.25 
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where g a is the spin-isospin degeneracy factor. The full shear viscosity and entropy density 
is a sum of the partial results, r\ = J^aVa and s = J2 a s a . We require that r] a /s a is equal 
to r]/s for each species, and thus the relaxation time parameter C® is, in principle, different 
for each species. However, for a classical gas the two integrals in square brackets are equal 
upon integrating by parts, and thus C® = r] a /s a = r\/s is independent of the mass and 
species label. For fermi-dirac and bose statistics these integrals are very nearly equal (to 
4% accuracy) for all values of the mass, and Ci is approximately equal to rj/s for all species 
independent of mass and statistics. 

Now that the relaxation time parameter C" is fixed for each species, the corresponding 
second order Sf for each species is found by appending a species label, n p — > n a v and C p — > 0°, 
to previous results. For a multi-component gas with a quadratic ansatz we find 



a 

V r n =Y,(C a p ) 



2g a 



p 
L(d-l)(d+l)(d + 3)T6y p ' l P E~ 



(d-l)(d + l)T*J p 



(3.8a) 
(3.8b) 



For classical statistics Cp 



rj/s, and integrating the first integral in square brackets by parts 
yields a simple relation noted previously [4] 



Ai = TjT^ 



(for a = 0) . 



(3.9) 



The remaining thermodynamic integrals are most easily done numerically; summing over all 
hadronic species with mass less than 1.5 GeV we find 



Ai = rjT^ 



rj 



(e + V) 



8.9. 



(3.10) 



Thus, TTn/(j]/s) = 8.9 would seem to be the most consistent value for the 2nd order trans- 
port coefficients during the hydrodynamic evolution of the hadronic phase. However, this 
value for Tr n is somewhat too large to be used comfortably in the simulation [19]. Further, 



There have been several efforts to go beyond this extreme species independent approximation [7, 18] 
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Tt-x decreases as temperature increases, and Tt^/(t]/s) ~ 5 is good approximation in the 
QGP phase [4]. We have therefore taken Ai = r]T w = 5r] 2 /(e + V) throughout the evolution. 
This means that there is a small inconsistency between the second order 5f at freezeout, 
and the second order parameters used to simulate the bulk of the hydrodynamic evolution. 
Similar inconsistencies are found in all attempts to consistently couple hydrodynamic codes 
with hadronic cascades [20]. 

To summarize, in this section we have specified precisely the initial conditions, the equa- 
tion of state, the transport coefficients at first and second order, and the first and second 
order corrections to the distribution functions. We have used this setup to compute the 
linear response coefficients v n /e n for RHIC and LHC initial conditions for the first six har- 
monics. Our results are displayed in Fig. 2 and Fig. 3. We will discuss the physics of these 
curves in the next section. 



IV. DISCUSSION 

This work computed the second order viscous correction to the thermal distribution 
function, 5/(2), and used this result to estimate the effect of second order corrections on the 
harmonic spectrum. Our principle results are shown in Fig. 1 for a conformal equation of 
state, and Fig. 2, and Fig. 3 for a lattice based equation of state. First, examine the t>2 
curves for the conformal EOS shown in Fig. 1. The most important remark is that even for 
a conformal equation of state, where the expansion is most violent, the derivative expansion 
converges acceptably for p? ~ 1-5 GeV, i.e. the second order correction is small compared to 
the first order correction. Not surprisingly, when a linear ansatz is used for Sf (rather than 
the more popular quadratic ansatz) the convergence of the derivative expansion is improved 
at high pt- Typically in hydrodynamic simulations of heavy ion collisions, the strictly 
first order 6fn\ is replaced by the 8f(\) which incorporates some, but not all, second order 
terms 3 . Examining Fig. 1, and also Fig. 2 and Fig. 3, we see that, while the sign of the second 
order correction is correctly reproduced by this incomplete treatment, the magnitude of the 
correction is generally significantly underestimated, and the pr dependence of the second 
order correction is qualitatively wrong. 

Most of these observations remain true for the more realistic lattice equation of state 
shown in Fig. 2 and Fig. 3. Generally, second order corrections are quite small for the 
first three harmonics, V\ to V3, and become increasingly important as the harmonic number 
increases. Indeed, for w 6 at RHIC and the LHC, the second order viscous correction is of 
order one, and the hydrodynamic estimate can no longer be trusted. It is also instructive 
to note that the sign of the second order viscous correction for p? ^ 1.5 GeV is positive, 
i.e. second order corrections bring the t> n (pr) curves closer to the ideal results. Generally, 
when the first order correction, becomes so large as to make v n {pr) negative at first order, 
the second order correction conspires to keep v n positive. When constraining rj/s with 
hydrodynamic simulations, the second order corrections are most important for V4 and v$. 
Indeed, at RHIC these corrections are quite important for v 5 even in central collisions. 

For a practical perspective, using Sf^) i n a hydrodynamic simulation is not particularly 
more difficult than using Sfm, and Eq. 2.18 can be readily implemented in most hydro 



As discussed above, Sf(i) uses IT** 1 ' in place of —r\o^ v when calculating the first order correction. 
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FIG. 2. Differential v n /e n from various viscous hydro dynamic simulations at b = 7.45 fm for RHIC 
and LHC initial conditions, and a lattice equation of state (LEOS) with Tf = 150 MeV. Here the 
n p curve shows the flow from second order hydrodynamics without the viscous correction to the 
distribution function; Sfm and 5fn\ + 5/(2) show the flow with the viscous correction at first and 
second order respectively; and finally, the 5fZ\ curve uses —rja^ instead of n^ u in the first order 
viscous correction (see Eq. (2.17)). 
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FIG. 3. Differential v n /e n from various viscous hydro dynamic simulations at b = 7.45 fm for RHIC 
and LHC initial conditions, and a lattice equation of state (LEOS) with Tf = 150 MeV. Here the 
n p curve shows the flow from second order hydrodynamics without the viscous correction to the 
distribution function; Sfm and 5fn\ + 5/(2) show the flow with the viscous correction at first and 
second order respectively; and finally, the 5fZ\ curve uses —rja^ instead of n^ u in the first order 
viscous correction (see Eq. (2.17)). 
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codes. The functional form of <5/( 2 ) and its magnitude is about as well constrained as 5/(i), 
and consistency with the second order hydrodynamic evolution would seem to mandate its 
use. At very least 5/(2) should be taken into consideration when estimating the uncertainty 
in the rj/s extracted from heavy ion collisions. Finally, when trying to use hydrodynamics in 
very small systems such as proton-nucleus collisions at RHIC and the LHC [21-27], second 
order corrections to Sf should be used in order to monitor the convergence of the gradient 
expansion. 
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Appendix A: Tensor decomposition of 5/(2) 

The goal of this appendix is to compute 5/(2) and to record how this results transforms 
under rotations in the local rest frame. Our starting point is Eqs. (2.18) which we rewrite 
in terms of irreducible tensors under rotations in the local rest frame. 

A systematic strategy decomposes all derivatives into temporal and spatial components 

d M = -u^D + V M , (Al) 

where the spatial component V M is orthogonal to u^. When differentiating a quantity that 
is already first order (i.e. P^d^Sf?^), we can use the lowest order conformal equations of 
motion to rewrite time derivatives in terms of spatial derivatives 

D\nT = -c 2 s V-u, (A2) 

D« M = -V A ,lnT, (A3) 

where c 2 s = (e + V)/(Tc v ) = l/(d— 1). Finally, the resulting tensors can be decomposed into 
symmetric, traceless, and spatial tensors as in Eq. (2.2), which transform irreducibly under 
rotations in the local rest frame. To illustrate the procedure, we record the decomposition 
of Da a p 

Da a p = D^A^A^) =A^A vfi Do^ + a^D(A m A vfi ) , (A4) 

= (Dcr aP ) - {<J^u a V„ In T + u p a v a V„ In T) . (A5) 

Similarly, the symmetrized spatial tensor {V^&ap},, m that arises when differentiating SfZ^ 
is decomposed as 

f 2 

{V M a Q/3 } sym = (V M a a/3 ) + | -j— -A^V^a} + u a (<r p fi a w ) 

+ u a -^-a 2 + 2u a (a p p n iMp )+2u a ^-V-u\ , (A6) 

U — 1 CL — 1 ) sym 



where we have used 



V m m p = -a MP + Q w + -^-V ■ u . (A7) 



2 



16 



Finally, we note that 

(9 A 7T A ) = A, X2 d Xl n x ^ = - V [(d - 2) <a MA V A In T> + <V A a A )] , 



(A8) 



where we have used the first order expression, tt^ = —rja 111 ', the conformal temperature 
dependence of 77 oc T^" 1 , and the lowest order equations of motion. 

With this automated set of steps, we start with Eq. (2.18) and place Sf7 2 \ into its canonical 
form 

pill ppu-i p/t 3 p^A p/l! p/l 2 pm 

&f(2) =Xi P ^ (^« ff »w) + X2 P ^5 [(V m a^ 3 ) - 3 M1M2 V M3 InT)] 



J>6 



y5 






4p 



-2 



d + 3 

pil2piH 

pi 



X2pE p 



\ p^piii 



pi 



!^c 



:A> 



( Dcr Wi) + 



a 



M2M1 
d-1 



V-w-2 



^aAia) 



+ 6*W [- <VaO - (d - 2) <a M2A V A lnT>] + ^ 



T 3 



pi 



(A9) 



where the functions Xop? Xip> X2 P and £ 3p and £ 4p are recorded in the text, Eq. (2.31). 
In this form it is easy to integrate over the phase space to determine the viscous stress 



-aw _ -aw 1 -Aw 



pupv 
P° 



(<% + <%) 



(A10) 



where 7T^ and 71"^ are given by static form of the constituent relation Eq. (2.5). Rotational 
invariance in the rest frame reduces these tensor integrals, e.g. 



pill pH2 pfl3 ptli 
XOp wj \^M3A*3/ 



p 

(d-l)(d+l).L X0P E r , 



(O^ 2 ) 



(All) 



yielding the equations for the transport coefficients written in the text, Eqs. (2.23), (2.28), 
and (2.29). In addition, we see that independent of the collision integral one finds the kinetic 
theory expectations identified in Ref. [4] 



A, 



-2777V , 



and 



A, = . 



(A12) 



Finally, in presenting these results in the text, and in implementing the results in a 
realistic hydrodynamic simulation, we have used the dynamic form of second order hydro- 
dynamics, where ii tJjU is treated as a dynamic variable. This choice amounts to using —tt^u/t] 
in place of <j^ v . In 8f? 2 ) this reparameterization yields the replacements: 



1 



(%»> + tttV • u - 2 «n wA ) ->. —K*» + ww\ 

(X J- ' I I TT 



1 



ill 



y^, 2 x) , (A13) 



(V w ^ w ) - 3 (<W 2 V M3 InT) -> - [{d + 2) (7T WM2 V M3 InT) - (V^vr^)] , (A14) 
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and Eq. (A8). In addition, when replacing a^ u with —7^/77 in the first order result, the 
difference -{j^u + Wv) must be appended to the second order result - see Eq. (2.19). The 
full result for 5 fa) and 5f{2) is given in Eqs. (2.20) and (2.30) respectively. 
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